#ifndef COMMON_CUDA_H
#define COMMON_CUDA_H

#include "common.h"

//calclulate the force between two bodies
inline __host__ __device__ float3 calcGravity(float4* pos_a, float4* pos_b)
{
	float3 r = make_float3(*pos_a) - make_float3(*pos_b);	// 3 flops
	float dis_2 = dot(r, r) + EPS;								// 5 flops
	float dis_6 = dis_2*dis_2*dis_2;						// 2 flops
	if (dis_6 == 0.0f)
		return make_float3(0.0f);
	else
		return (-G*(pos_a->w*pos_b->w)*rsqrtf(dis_6))*r;	// 7 flops
}

//calclulate the energy between two bodies
inline __host__ __device__ float calcGravityPot(float4* pos_a, float4* pos_b)
{
	float3 r = make_float3(*pos_a) - make_float3(*pos_b);	// 3 flops
	float dis_2 = dot(r, r);								// 5 flops
	if (dis_2 == 0.0f)
		return 0.0f;
	else
		return 0.5f*(-G*(pos_a->w*pos_b->w)*rsqrtf(dis_2));	// 
}

//calclulate the acceleration of a to b
inline __host__ __device__ float3 calcGravityAccel(float4* pos_a, float4* pos_b)
{
	float3 r = make_float3(*pos_a) - make_float3(*pos_b);	// 3 flops
/*
	float dis_2 = dot(r, r);								// 5 flops
	float dis_6 = dis_2*dis_2*dis_2;						// 2 flops

	if (dis_6 == 0.0f)
		return make_float3(0.0f);
	else
		return (-G*(pos_b->w)*rsqrtf(dis_6))*r;					// 6 flops
*/

	float inv_dis = rsqrtf(dot(r, r) + EPS);
	float inv_dis_3 = inv_dis*inv_dis*inv_dis;

	return (-G*(pos_b->w)*inv_dis_3)*r;					// 5 flops
}


//update the positions
__global__ void updatePos(float4* pos, float4* vel, int n, float dt);
/*
{
    int x = blockIdx.x*blockDim.x + threadIdx.x;

    if (x < n)
        pos[x] += dt*vel[x]; // 3 flops
}
*/

#define SUM_BLOCK_SIZE 128

template <class T>
__global__ void foldSum(T* f, int n, int pivot)
{
	int x = blockIdx.x*SUM_BLOCK_SIZE + threadIdx.x;

	if (x < n && x >= pivot)
		f[x - pivot] += f[x];
}

template <class T>
void parallelSum(T* a, int size)
{
	dim3 dimb(SUM_BLOCK_SIZE, 1);
	dim3 dimg((size-1)/dimb.x + 1, 1);

	int pivot = msb32(size-1);
	int curr_size = size;

	int i;
	for (i = pivot; i > 0; i >>= 1)
	{
		foldSum<<<dimg, dimb>>>(a, curr_size, i);
		curr_size = i;
	}
}

#endif
